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When relaxation towards an equilibrium or steady state is exponential at large times, one usually con¬ 
siders that the associated relaxation time t, i.e., the inverse of that decay rate, is the longest characteristic 
time in the system. However that need not be true, and in particular other times such as the lifetime of an 
infinitesimal perturbation can be much longer. In the present work we demonstrate that this paradoxical 
property can arise even in quite simple systems such as a chain of reactions obeying mass action kinetics. By 
mathematical analysis of simple reaction networks, we pin-point the reason why the standard relaxation 
time does not provide relevant information on the potentially long transient times of typical infinitesi¬ 
mal perturbations. Overall, we consider four characteristic times and study their behavior in both simple 
chains and in more complex reaction networks taken from the publicly available database “Biomodels”. In 
all these systems involving mass action rates, Michaelis-Menten reversible kinetics, or phenomenological 
laws for reaction rates, we find that the characteristic times corresponding to lifetimes of tracers and of 
concentration perturbations can be much longer than t. 

PACS numbers: 


I. INTRODUCTION 

Networks have been used to model systems involving 
large numbers of components, agents, or species UTJ. Of 
particular interest are the effects arising in such systems 
either because of out-of-equilibrium dynamics or through 
equilibrium phase transitions. Collective effects are gener¬ 
ally associated with slow dynamics, i.e., characteristic times 
that are much larger than the microscopic times associated 
with elementary processes. In the present work our focus 
is on the emergence of large characteristic times in reac¬ 
tion networks close to their steady state. There are many 
ways to define a characteristic time in a dynamical system. 
The simplest is via the asymptotic relaxation towards the 
steady state [HO, relaxation which often will be exponen¬ 
tial. If so, the amplitude of the perturbation or “distance” 
to the steady state will decay as exp(— t/T) at very long 
times, from which one then defines t to be the relaxation 
time. Although in familiar situations t is the longest char¬ 
acteristic time, our goal here is to investigate cases where 
much larger times can arise. Our study focuses on reaction 
networks for specificity, but our framework is generally ap¬ 
plicable to any system. 

Reaction networks involve species that can transform 
one into another. If the species are molecular, one can 
get insights into the dynamics of the system by introduc¬ 
ing an isotopic tracer and by following in time its incor¬ 
poration into the different molecular species [j4|]. Assume 
that the reaction network is in contact with outside reser¬ 
voirs, and let t t be the time the tracer takes to exit the 
system. Surprisingly, the mean of t t , corresponding to the 
lifetime [ElO] of the tracer (and sometimes called the mean 
residence time of the tracer), can be much greater than t. 
The object of our work is to understand such a possibil¬ 
ity, pointing in particular to the danger of assuming that 
t is the main and longest characteristic time in these sys¬ 
tems. For pedagogical reasons, we will begin by treating 
one-dimensional networks because an in-depth analytical 
treatment is feasible there, from which one can easily un¬ 


derstand the influence of network size. We will then study 
more general systems using reaction networks published 
by other authors. In all cases, we compare the behaviors of 
four characteristic times in these systems, investigating the 
causes that can render them non informative or make their 
ratios diverge. 


II. MODELS AND METHODS 

A. Networks, molecular species and associated reactions 

A metabolic network consists of a set of reactions and 
associated metabolites. It is convenient to represent such 
a network as a graph where the nodes are associated with 
metabolites; these are linked together by edges when there 
is a reaction that includes them as substrate and product. 
Such edges may be uni or bi-directional, accounting for 
the reversibility of the associated reaction. Let there be N 
metabolites M t (i = 1, ...,1V) and define C, as the concen¬ 
tration of M;. We are interested in the dynamics of the C h 
i.e., how these quantities change with time and in the cor¬ 
responding fluxes through the different reactions. Specif¬ 
ically, we shall study the dynamics close to the system’s 
steady state and we shall probe the associated character¬ 
istic times. To facilitate the mathematical understanding 
of these times, we shall first focus on a particular kind of 
network consisting of a linear chain of reactions. In that sit¬ 
uation, we order the metabolites from 0 to N +1 where the 
metabolite M, is the product of reaction R t whose substrate 
is metabolite M l _ 1 : 

Vi Vo V W V/U4.-1 

M 0 +1* Mj ... M n M n+1 (1) 

The metabolites M 0 and M w+1 reside in infinite reser¬ 
voirs at the two extremities of the chain so their concen¬ 
trations are constant. By convention, the forward direction 
in such a chain goes from M 0 to M N+1 . Once understood 
the characteristic times in this system, we shall use the in- 
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sight thereby gained to probe the situation in more realistic 
metabolic networks having branches and loops. 

Reactions transform metabolites into other metabolites 
but it is necessary still to specify the actual kinetics. When 
a reaction happens spontaneously, without the need for a 
catalyst, it can be modeled by a mass action rate law (MA) 
where the flux is given by 

vf iA = a i C i _ 1 -b i C i . (2) 

To be specific, one can consider using the usual conven¬ 
tion whereby concentrations are measured in Moles per 
liter and fluxes in Moles per liter per second. The param¬ 
eter a t (resp. b,) is then the probability per second that 
a molecule of metabolite (resp. spontaneously 
transforms into a molecule of metabolite M i (resp. M,^). 
Note that Eq.[2]gives the total flux which is the forward flux 
minus the backward flux. 

In practice, one is often interested in catalyzed reactions 
where the spontaneous rates are terribly low. For instance, 
in biochemistry, most reactions are catalyzed by enzymes; 
the catalysis allows for rates that can be enhanced by a 
factor of 10 10 or more. For any such enzymatic reaction, 
the rate may be limited by the amount of enzyme and is 
no longer directly proportional to metabolite concentra¬ 
tion. Generally, the relation between substrate concentra¬ 
tion and reaction rate grows linearly at low concentrations 
and then saturates at high concentrations of substrate. The 
reaction kinetics in this situation are typically modeled by 
the so called reversible Michaelis-Menten-Henri (MMiT) 
law ED- In the case of a reaction involving one substrate 
and one product, the flux is given by 
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Here, a, is the maximum rate in the forward direction, 
reached when the substrate is in large excess and the prod¬ 
uct is absent. Similarly, /3 ; is the maximum rate in the 
backward direction. The maximum forward rate is pro¬ 
portional to the enzyme concentration and is often decom¬ 
posed as a = k cat E with E being the enzyme concentration 
and k cat the maximum number of reactions catalyzed by 
one molecule of enzyme per unit of time. K l s and K l p , called 
the Michaelis constants respectively for substrate and prod¬ 
uct, are characteristic concentrations which set the scale 
for when the reaction becomes saturated in substrate or in 
product. For a MMH reaction in the absence of the prod¬ 
uct, K s is the concentration for which the rate is at half of 
its maximum value. 


equal, corresponding to tuning the concentrations so that 
their ratio is the equilibrium one. Outside of that special 
case, the system will be out of equilibrium and concentra¬ 
tions will change in time until a steady state is reached 
which necessarily will have non zero fluxes. This steady 
state is generally unique if there are no regulatory pro¬ 
cesses but for our study to be completely general, we will 
not assume uniqueness of the steady state, we shall simply 
consider a stable steady state and investigate its character¬ 
istic times. 

We have followed two approaches for determining 
steady states (leading to identical results): 

1. solve the steady state equations dCjdt = 0 which 
we do numerically using root finding (routine “find- 
root” in Python). For any given boundary conditions, 
i.e., concentrations C 0 and C N+1 , this leads to a list 

of steady-state concentrations C ss . It is necessary to 
check that the resulting steady state is linearly sta¬ 
ble. This check can be performed using the linearized 

equations about the steady state. If 5C is the (in¬ 
finitesimal) difference between the actual concentra¬ 
tions and those in the steady state, one has 


d5C 

~dt~ 


= J (c) 5C 


(4) 



A; if j = i — 1 
_ — QVkl + B i) if j = i 
B i+ 1 if j = i + 1 
0 otherwise 


(5) 


where the A, and B, are related to the terms entering 
Eq.[2]for mass action and Eq.[3]for Michaelis Menten 
Henri as specified in Table]!] is the NxN Jacobian 
matrix with indices i and j going from 1 to N; the 
superscript c refers to the fact that it describes the 
(linearized) dynamics of (perturbed) concentrations. 
The steady state is stable if all the eigenvalues of the 
Jacobian have negative real part. 

2. follow the concentrations using the dynamical equa¬ 
tions (the system of ordinary differential equations 
specified by the kinetic laws) and extract the long 
time limit of the concentrations. This requires extrap¬ 
olation, but generically takes one to a stable steady 
state. 


C. Defining four characteristic times 


B. Determining steady states 

When a physical system is not driven by outside forces, 
it goes to its equilibrium state where all net reaction fluxes 
are 0. In the context of our one dimensional model, that 
can only arise if the free energies of the two reservoirs are 


• The first characteristic time is the relaxation time de¬ 
fined as — l/A^ where h[ c> is the real part of the 
leading eigenvalue of J (c) having the largest real part. 
Because this time is defined via the linearized dynam¬ 
ics for the concentrations about the steady state, we 
shall refer to it as t c . 
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• The second characteristic time is the previously men¬ 
tioned tracer lifetime (or mean residence time), 
which we denote by T t . The motivation for intro¬ 
ducing this quantity comes from tracer experiments 
in chemical networks where isotopic labels are used 
to follow atoms as reactions progress. Instead of in¬ 
troducing a perturbation to concentrations, this ap¬ 
proach labels (e.g. via an NMR pulse) atoms of one 
metabolite M k at t = 0 without changing any con¬ 
centrations. In practice this labeling affects only a 
fraction of the molecules. The effect of this labeling 
is to leave the fluxes unperturbed as well. The system 
stays in its steady state, it is just that some of these 
concentrations become labeled. Note that when one 
labeled metabolite is tranformed into another, the la¬ 
beling follows because the labeled atoms. 

Let us study the time evolution of the concentrations 
of these tracers C t = {C t l ,C t2 ,---,C tJV } (the sub¬ 
script t is for tracer ). As previously introduced, let 

C ss = {q s ,q s ,...,C“} be the steady state concen¬ 
trations. Consider the reaction R, and let <p f ( be its 
forward flux and cp b i its backward flux in the steady 
state. Then the labeled concentration C t i will include 
an incoming term given by the rescaled forward flux 
c/jf ; C t i _ 1 /C s i s _ 1 because all metabolite molecules (la¬ 
beled or not) have an equal probability of participat¬ 
ing in the reaction R t . As a result, the dynamics of 
the tracer concentrations is 

dC t > 

- r £ = J (t) C t (6) 

dt 


• The previous definition of lifetime of a tracer can be 
generalized to the lifetime of any quantity and in par¬ 
ticular to a perturbation to steady-state concentra¬ 
tions. Suppose one introduces at t = 0 an infinitesi¬ 
mal perturbation in the concentrations, 5C(0). Then 

according to Eq. 0 <5C(t) = exp(tJ (c) )5C(0). In di¬ 
rect analogy with Eq. [8j the lifetime of that perturba¬ 
tion is 


, S™\5CW\dt 

|5C(0)| 


(9) 


providing a third characteristic time of our system, 
referred to as the lifetime of a concentration per¬ 
turbation. To be completely general, both here and 
for the tracer lifetimes, the vectors of concentrations 
should be actually the deviations of their values from 
their long time limit. Indeed, if there were no reser¬ 
voir and thus no exit possible of the atoms, the long 
time limit of the perturbation or tracer concentration 
would not be 0. 


• Our fourth and last characteristic time is z t , defined 
as — 1 / A ^ where A^ is here the real part of the lead¬ 
ing eigenvalue of J {t) . It corresponds thus to the usual 
relaxation time but for the tracer molecules rather 
than for the metabolite concentrations, thus the sub¬ 
script t. 


III. BEHAVIOR OF CHARACTERISTIC TIMES IN THE 
ONE-DIMENSIONAL NETWORK 



>f/Cf if j = i -1 

« -(<!>{/qq + tljqq) if j = i 
4> b jq s if j = i +1 

1 0 otherwise 


(7) 


Note that these linear dynamics are exact even if C t; 
is not infinitesimal. In general, the matrix has no 
reason to be identical to J w . By exponentiating, one 
has the expression for the labeled concentrations at 
all times: C t (t) = exp(tJ M )c t (0). The lifetime T t of 
the tracer is then obtained as the average over time 
of the survival probability: 


fo\qM\dt 

|Q(0)| 


( 8 ) 


In this equation, |C t (t)| is the norm of the vector. 
For our study, we use the L 1 norm because it makes 
more sense for an atomic tracer which is conserved. 
Note also that Tf in Eq. [8] is the direct analog of the 
mean lifetime of a decaying positive scalar quantity; 
the norm allows one to extend the notion to a vector 
in a straightforward manner. 


As can be seen from the four characteristic times defined 
in the previous section, we distinguish two properties of 
a metabolic system: (i) the dynamics of an infinitesimal 
perturbation in the concentration of metabolites and (ii) 
the spreading and drift of tracers. Each of these proper¬ 
ties can be considered when reaction kinetics are given by 
MA or MMH rate laws. In each case one can define both 
the standard relaxation time based on the asymptotic decay 
rate and a lifetime which measures the characteristic time 
needed for the system to return close to its steady state. 
In the case of a chain of reactions with the same kinetic 
parameters, the homogeneity allows us to obtain results 
analytically. For instance in the case of MA, the linearized 
dynamics (J (c! and J (t) ) are independent of the steady state 
chosen (that is the concentrations of M 0 and M N+1 do not 
enter) and the matrices are sufficiently simple for one to 
obtain in closed form the eigenvectors and eigenvalues. In 
the case of a MMH framework, when one performs the lin¬ 
earization about the steady state, the resulting system is ho¬ 
mogeneous only if the steady state itself is homogeneous, 
which requires that all the metabolites have the same con¬ 
centrations. When this is the case, the steady state is again 
obtained in closed form. Furthermore, the eigenvectors 
and eigenvalues can be derived analytically, which gives us 
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Position 


Figure 1: Decrease with time of |5C| (or | C t |) for an initial 
t = 0 excess concentration (or labeling) localized at a site 
in the middle of the chain of reactions. The y axis is on a 
log scale so that one can see the asymptotic exponential 
decay as a straight line of slope — 1/t: t 20 = 4.92, 
t 50 = 5.65, and t 100 = 5.78. All N mass action reactions 
have a = 2 and b = 1. Shown are cases with N = 20,50 
and N = 100. 


then the formulas for t c and r t . Unfortunately the study 
of the lifetimes T c and T t requires resorting to numerical 
methods, but these are relatively straightforward as they 
reduce to calculating exponentials of the matrices J fcj and 
J- 4 - 1 and performing the integrations in Eq.j9jand[8] For the 

initial perturbation, for simplicity we take 5C(0) and C t (0) 
to vanish everywhere except on the site at the center of the 
chain where the value is set to 1. For an even number of 
sites, there is no such center so we average over the two 
most central sites. 


A. Long transient times drive the gap between lifetimes and 
relaxation times 

The integral in Eq. [^depends on C £ (t) = exp(tJ (t) )C t (0) 
which can be written using spectral decomposition as a sum 
of N terms, each term being associated with one eigenmode 
and having the time dependence exp(tA £ c) ) where is 

the associated eigenvalue. When N = 1, C t (t) is a con¬ 
stant times a single decaying exponential. Plugging into 
Eq. [8] then reveals that T t = z t . The paradox whereby T t 
can be much larger than r t arises only when NA> 1. It 
is true that each of the N terms contributing to the spec¬ 
tral decomposition of C t (t) decays in magnitude at least 
as fast as exp(— t/r) but that does not mean that the sum 
of these terms has that behavior on time scales comparable 
to t. Indeed, the terms are not all of the same sign, and 
their cancellations can lead to long transients before the 
asymptotic behavior (the exponential decay) prevails. To 
illustrate this, we show in Fig. [l] the L 1 norm of C t (t) as a 
function of t in our toy model consisting of a chain with a’s 
and V s identical across MA reactions. At large times, one 


Figure 2: Leading eigenmode profile for the 20 last 
metabolites of the chain. Parameters: a = 2, b = 1, and 
N = 20,50,100. 


sees the exponential decay (a straight line on this semi-log 
plot) but this asymptotic behavior may set in at times only 
much longer than t itself. The cancellation at short times 
just mentioned is particularly striking: the curve is very flat 
for a very long time before it begins to decrease. That wait¬ 
ing time contributes to the large difference between T t and 
T t and is associated with the transient time one must wait 
for tracer molecules to exit the system. 


B. Dependence of the characteristic times on N 


Assuming the reactions to all have the same parameters 
and that the steady state is also homogeneous (cf. previ¬ 
ous remarks), the relaxation time (be it t c or T t ) can be 
obtained by using the translation invariance of J {c) and J (LJ . 
Each eigenvector is a product of a sine and an exponential. 
The formula for the eigenvalues leads to 

1 

T = - - - 

A + B — 2-J AB cos 



where the quantities A and B are the forward and backward 
probability of transition per unit of time in the equations 
linearized about the steady state, entering in J (c) for t c and 
in J (t) for T t . They depend on whether one considers MA or 
MMH reaction kinetics and whether one considers a con¬ 
centration perturbation or a tracer, the different cases being 
enumerated in Tabled 

The ts in the four cases are given by a standardized for¬ 
mula (Eq. 101, it is just that the proper A and B coeffi¬ 
cients must be used. Note that for MA kinetics, J {c:) = J® so 
t c = T t . Furthermore, in both MA and MMH frameworks, 
when the relative difference between A and B is small, the 
ts exhibit two different regimes, one for small chains and 
one for long chains. For a short chain, N -C N cross = 
the characteristic times grow quadratically with the num¬ 
ber of metabolites in the chain, a feature characteristic 
of diffusing systems for the simple reason that if A = B, 
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Figure 3: Relaxation times {left) and lifetimes (right) for chains between 2 and 100 metabolites long for a perturbation of 
concentrations and a tracer using the mass action or the Michaelis-Menten-Henri framework. Parameters: a = 2 and 


b = 1, K s = K P = 2 and a = aK s and ft = bK P so that the three conditions are comparable. The large N relaxation times 
are respectively T^f im = t ^f im = 5.83, = 11.66, = 47.66. The transition sizes between a quadratic and 


constant or linear behavior are 

c, cross t, cross 


i-i -ktMMH _ y lyMMH _ i y 

't, cross ^c, cross ' 


Parameter 

MA-c MA- 

t MMH- 

cMMH-t 

A 

a a 

a-F 

a 



K S S 

K S S 

B 

b b 

P+F 

_ 2 _ 



K P S 

K P S 


Table I: Value of the A and B parameters for the four 
situations considered. F and S = (1 + c ss /K s + c ss /K P ) are 
respectively the flux and the saturation factor at steady 
state in the network for the reactions, the system being by 
hypothesis homogeneous. The “c” (respectively the “t”) 
appended to MA and MMH denotes that it is the 
perturbed concentrations (respectively the tracer 
concentrations) that are concerned. 


the dynamics is purely diffusive. When N is much above 
this crossover value, t c and r t become independent of the 
chain length as can be seen directly by setting to 1 the co¬ 
sine in Eq. [l0} 

Note that the crossover size N cross diverges as the inverse 
of A — B. Furthermore, in the context of MMH reaction ki¬ 
netics, this crossover occurs for larger chain lengths when 
considering the dynamics of a concentration perturbation 
than when considering tracers because the saturation has 
the effect of reducing the difference between A and B. To 
illustrate these effects, we display in Fig. [3] the relaxation 
times as a function of the chain length N for particular val¬ 
ues of the kinetic parameters. As for MA, t c and T t do not 
increase asymptotically with N, the characteristic times be¬ 
come independent of the system size. To understand how 
this occurs, let us examine the leading eigenvector. Its en¬ 
tries depend exponentially on the index i of the node and so 
its profile is biased towards the largest indices. If the eigen¬ 
vector with the largest eigenvalue becomes dominant, the 
major part of the deviation from the steady state is located 
on a few metabolites (about N cross ) at the end of the net¬ 
work. As illustrated in Fig.|2j if one increases the number of 


metabolites, that eigenmode just gets shifted to stay at the 
same position when measured from the end of the chain. 
As a consequence, increasing N does not affect the corre¬ 
sponding eigenvalue which determines t. Thus r c and r t 
become independent of N at large N. 

For the T c and T t lifetimes, we did not derive a closed 
form expression but one can still distinguish between two 
regimes. If A — B is small the behavior for small N is 
diffusion-like so T c and T t increases quadratically with N. 
In contrast, for long chains, if A / B, one has a regime 
where T c and T t grow linearly with N. Similar arguments 
as for the relaxation times t can be invoked to explain these 
two regimes. In small networks, the diffusion to the two 
sides of the chain dominates over the mean drift toward 
one end of the chain. In large networks, assuming A > B, 
most of the transient time dominating T c and T t is dedi¬ 
cated to the transport of the molecules to the N + 1 end, 
therefore that transient time is roughly equal to N divided 
by the drift velocity (which is proportional to (A — B)). We 
illustrate these behaviors in Fig. [3] where one sees again 
that the various cases behave similarly with the network 
length. (We already noted that for MA kinetics, J {c:) = J®; 
as a consequence one has T c = T t there, just as one has 
T c = T t .) 


C. Effect of the saturation on the characteristic times 

The major differences between MA and MMH come 
from the effect of the saturation. In the case of the MA rate 
laws, there is no saturation while saturation effects can be 
important in MMH kinetics. This difference can lead to 
much larger characteristic time scales in MMH than in MA 
whenever the concentrations are larger than K s or K P . Fur¬ 
thermore, for highly saturated enzymes, the characteristic 
times can be very different depending on whether one ob- 
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Figure 4: Relaxation times and lifetimes as a function of 
the saturation. Parameters: a = 4, [3 = 2, K P = 2, N = 30. 
To vary the saturations, the parameter K s is changed over 
a range going from 1 to 10 -4 . 


serves a tracer or a perturbation of concentration. Consider 
a reaction that is near saturation. Introducing a perturba¬ 
tion in the substrate will not change much the flux of that 
reaction and as a result it will take a long time to dissi¬ 
pate the perturbation away. On the other hand a tracer is 
essentially unaffected by saturation effects. Indeed, it is 
not because the reaction is saturated that the tracers can¬ 
not participate in the reactions. In effect, the tracers freely 
pass reactions that are saturated. The main consequence of 
this phenomenon is that in MMH r c can be much larger 
than r t (and T c can be much larger than T t ). 


To investigate quantitatively this phenomenon of partic¬ 
ular relevance when interpreting kinetic properties from 
tracer measurements, let us increase saturation effects by 
reducing K s . K P could also have been reduced, but when 
doing so, the flux in the network may reverse which un¬ 
necessarily complicates the analysis. Using the parameters 
of Table [I] in Eq. 10 for small values of K s gives the fol¬ 
lowing analytical values for the dominant terms of the two 
relaxation times associated with a tracer (T t ) and with a 
concentration perturbation (t c ): 



IV BEHAVIOR OF CHARACTERISTIC TIMES IN MORE 
GENERAL METABOLIC NETWORKS 


A. Effects of disorder in the one dimensional chain 

In the disordered (i.e., heterogeneous) case we now con¬ 
sider, the rates “a” and “b” for the different reactions are 
taken to be independent random variables. Because every 
rate is a positive variable, we draw it from a lognormal 
distribution, i.e., the natural logarithm of a rate r ; is dis¬ 
tributed according to a Gaussian of mean p and standard 
deviation a. Consequently, the mean of r { is p = exp(,u + 
a 2 / 2) and its variance is d 2 = (exp(cr 2 ) — 1) exp(2/i + a 2 ). 
We impose p to be equal to the value of the rate in the ho¬ 
mogeneous case. An appealing feature of that way of intro¬ 
ducing disorder is that the mean drift velocity of a marked 
molecule in Mass Action remains unchanged, being equal 
to its disorder average, (a t — b ; ). We are then left with the 
parameter d which can go from 0 to oo and quantifies the 
intensity of the disorder. In practice, we use the same coef¬ 
ficient of variation for the “on” and the “off” reaction rates, 
corresponding to a single measure of intensity of disorder: 
CV = a a /a = a b /b. 

For weak disorder, one expects little change in the values 
of the characteristic times (t c ,t t ,T c ,T t ) compared to the 
homogeneous case. However, as disorder (CV) increases, 
the characteristic times typically increase. To identify the 
typical behavior, we have determined these characteristic 
times for 10,000 realizations of the disorder and calculated 
the median times. We illustrate the associated results in 
Fig. [5]for t c and T t in the case of Mass Action where those 
two quantities are equal. Increase is relatively mild (cf. the 
scales) at low CV butis more marked when CV is larger 
than 30%. 

Consider now the effects of disorder on the lifetimes. In 
Mass Action, T c = T t , even in the presence of disorder. We 
display in Fig. [5] the dependence of these quantities on N 
for several values of CV and see that disorder has little ef¬ 
fect as long as CV is small. This can be justified by noticing 
that the drift velocity of a molecule at site i is a ; — b i _ 1 and 
its ensemble average (as in an annealed approximation) is 
the same as without disorder, namely a —b. At large disor¬ 
der this argument fails because the quenched and annealed 
averages are very different. An extreme case can be seen 
from the fact that a large value of “a” at one site cannot 
compensate a small value at another site. At large CV, one 
sees clear effects of disorder. The reason should be clear: T c 
and T t are sensitive to unfavorable reactions (for instance 
where a is small) throughout the whole chain of reactions. 


We see from these equations that r t is independent of the 
saturation while r c behaves linearly with 1/K S . Note that 
the saturation S = 1 + c ss /K s + c ss /K P scales in the same 
way for small K s . In Fig. [4] we show the dependence of the 
ts and the Ts on the saturation S for both a tracer and a 
concentration perturbation, assuming MMH rate laws. Not 
surprisingly, T c is strongly affected by S, just as t c is. 


B. Networks with branches and loops 

Although quite a few biosynthetic pathways include suc¬ 
cessive steps forming a chain of enzymatic reactions, the 
one dimensional systems considered so far remain toy mod¬ 
els because in all known organisms, large scale biochemi¬ 
cal metabolic networks have numerous branches and loops. 
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Figure 5: Relaxation (left) and transiting/it) times as a function of N for several intensities of disorder in the reaction rates 
as measured by their coefficient of variation CV. Main figure: T = T c = T t in the case of Mass Action kinetics. The error 
bars show the 68.2% confidence interval, value motivated by taking one standard deviation on both sides of the mean of 
a Gaussian distribution. Parameters: a = 2, b = 1, CV = 30% and 70%. 


It is thus necessary to consider how characteristic time 
scales might be affected by such structures. Rather than 
produce artificial networks including those features, it is 
more relevant to study directly the various kinetic models 
of metabolism that have been proposed in the literature. 
The repository “Biomodels” MM provides the gold stan¬ 
dards for such models both because the models must past 
tests to be deposited and because their availability ensures 
that they can be compared to state of the art. Focusing fur¬ 
ther on those models that have been manually curated, we 
are left with only a handful of cases. The reason is that 
measuring kinetic constants of enzymes is a very difficult 
task so almost always when building a kinetic model the 
modeler has to use indirect methods to overcome the prob¬ 
lem of dealing with many unknown parameters. We stud¬ 
ied four of these models, published respectively in 1110441311 . 

For each of those four kinetic models, we first down¬ 
loaded its SBML specification Q9] from the repository and 
exported the ordinary differential equations into Python 
code that can be processed. Once in our format, we de¬ 
termined the steady state of the network of reactions and 
we then computed the matrices J (c:) and J®. The associ¬ 
ated leading eigenvectors and eigenvalues were obtained 
using the inverse power method, thereby providing the val¬ 
ues of t c and T t . Furthermore numerical integration was 
used to compute T c and T t according to Eqs. [8jmd[9] The 
initial perturbation was taken to be localized at the first 
metabolite produced from the compound entering the net¬ 
work from the outside reservoir. 

In Table [IT] we provide the values of the four character¬ 
istic times for each of the Biomodels studied. The first 
model IllOll contains the reactions for glycolysis in S. cere- 
visiae (baker’s yeast). It has 17 reactions, mostly of the 
reversible MMH type, and there are 14 internal metabo¬ 
lites. Glucose is an external metabolite which enters the 
metabolism and then gets transformed. A total of 3 com¬ 
pounds can be excreted, all irreversibly. The characteris¬ 


tic times of this model are short, from a few seconds to a 
few minutes. Further inspection shows that the ordering 
of these four values follows the same pattern as in our one 
dimensional toy model, namely 

r t < r c < T t < T c (12) 

This can be justified as follows. First, r t > r c and T t > T c 
because a labeled atom is not subject to Michelis-Menten 
saturation effects. The saturation of flux in a reaction may 
prevent a concentration fluctuation from being evacuated 
but it will not prevent labeled atoms from going through 
(participating to the flux). Furthermore, in our toy model, 
the ts are relatively insensitive to processes inside the net¬ 
work, they depend mainly on reactions close to the ex¬ 
creted metabolites, while the T s depend on drift through¬ 
out the whole network and thus should be larger than the 


The other models follow quite closely this same pattern 
(cf. Table Model 2 contains the reactions for the gly¬ 
colysis and the pentose phosphate pathway in E. coli HlHI . 
It has 48 reactions and 17 internal metabolites, but we 
needed to remove the model’s explicit time dependence 
to obtain a meaningful steady state. The main difference 
with the model 1 is the modelled organism and the glu¬ 
cose steady state uptake rate (3.1 fimol.s -1 .L -1 compared 
to 1.5 mmol .s -1 ,L _1 ) but Eq. 12 is respected. Again model 


3 contains the glycolysis and the pentose phosphate path¬ 
way, but for a human cancer cell. It has 29 reactions and 
34 internal metabolites. The glucose uptake, expressed per 
gram of cell dryweight (0.17mmoZ.s _1 .gcdw _1 ), cannot be 
compared to the two previous uptakes but Eq.[l2]is mostly 
verified. 

Model 4 contains the reactions for the biosynthesis of 
purines in E. coli 111311 . It has a total of 29 reactions and 18 
internal metabolites. The main difference compared to the 
other three models is that the formalism uses kinetics that 
are neither MA nor MMH: the forward and backward rates 
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of the reactions are fractional powers of the concentrations 
of the metabolites. Such fractional powers are often used 
phenomenologically to parametrize allosteric or regulatory 
effects; they have the drawback that the flux may rise very 
steeply when starting with low concentrations; although 
this may be the case for some regulatory processes, it can 
lead to a situation where a concentration perturbation will 
be evacuated more efficiently than a labeled atom. Such a 
possibility seems to be realized in this model as in Table [II] 
one sees that t c < T t and T c <T t . 


time (s) 

Tc 

T t 

T c 

T t 

Model 1 m 

15. 

3.75 

339 

84.4 

Model 2 OH 

120 

95.2 

2834 

2210 

Model 3 [jT2| 

4.94 

0.16 

107 

3.53 

Model 4 G3 

4.34 10 1 2 3 4 5 

1.11 10 6 

9.35 10 6 

2.36 10 7 8 9 


Table II: Value of the characteristic times t c , T t , T c and T t 
in seconds for the four manually curated models 1110441311 
we have studied and that are available on the Biomodels 
repository [J8] • 


V DISCUSSION AND CONCLUSIONS 

The damping of a concentration fluctuation generally re¬ 
quires the perturbation to spread out but in our reaction 
network the drift plays a central role. The time scale for 
evacuating a perturbation is what we call its lifetime T (cf. 
Eqs.[8]and [ 9 ]), though in other contexts it can be referred 
to as the mean residence or transit time. In the absence of 
a drift, corresponding to a pure diffusive regime, the life¬ 
time T scales as the square of the diameter of the network, 
a scaling which also arises for the standard measure of re¬ 
turn to equilibrium via the relaxation time t. That is the 
situation one is most familiar with, and there t provides 
the longest characteristic time as it should. However, for 


typical reaction networks, one has both diffusion and drift. 
In particular, out of equilibrium systems will have fluxes, 
and such fluxes may drive labeled atoms out of the network 
just as drift does. In the presence of such drift, a perturba¬ 
tion’s lifetime T can scale as the diameter of the network 
divided by a characteristic drift velocity which is related to 
the presence of flux. Interestingly, in this out of equilib¬ 
rium situation, the relaxation time t is no longer informa¬ 
tive about the time scale of the (slow) process which evac¬ 
uates perturbations. In particular, in our toy model consist¬ 
ing of a homogeneous chain of reactions, t did not grow 
with the system size while T grew linearly. We showed 
analytically how that could be in that system, but the phe¬ 
nomenon is general. Indeed, in the presence of drift, the 
linearized dynamics can be decomposed into eigenvectors, 
but the leading eigenvector determining t tends to be con¬ 
centrated on the metabolites that can be excreted. As a re¬ 
sult, t is quite insensitive to the size of the network while 
T inevitably increases with network size since the evacu¬ 
ation of a perturbation requires it to cross the diameter of 
the network. These phenomena are most easily understood 
when the reactions obey mass action, but they arise also for 
Michaelis-Menten-Henri reaction laws. For this last case, 
the existence of a saturation of the flux with concentration 
of metabolites exacerbates the difference between T and 
t. Interestingly, the dynamics of labeled atoms that are of¬ 
ten used to investigate kinetic properties of networks are 
far less sensitive to these saturation effects. As a conse¬ 
quence, the use of isotopic labelings can lead to severely 
underestimate the longest characteristic time in these reac¬ 
tion networks. 
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